58. Multigrid and Multilevel Methods

Multigrid (MG) and Multilevel (ML) algorithms provide preconditioners with optimal condition numbers \(\kappa (C^{-1} A) = O(1)\), and optimal computational complexity \(O(N)\).

They can be seen as extension of the two level overlapping domain decomposition method to more levels.

• Ulrich Trottenberg, Cornelius W. Oosterlee, Anton Schuller: Multigrid, Academic Press, 2001

• Wolfgang Hackbusch: Multi-Grid Methods and Applications, Springer, 1985

In short, both methods are sub-space correction methods where the space splitting is defined by all basis functions together from a hierarchy of grids. Multilevel methods are additive Schwarz methods, while Multigrid methods are the multiplicative Schwarz methods. Both can be implemented efficiently by recursive algorithms.

We will present different theories for their analysis, one is based on sub-space decomposition using the ASM lemma, the other one uses the smoothing-and-approximation properties.

58.1. Multilevel preconditioner

We have a sequence of hierarchically refined meshes, which lead to a sequence of nested finite element spaces

\[ V_0 \subset V_1 \subset \ldots V_L \]

of dimension \(N_l = \operatorname{dim} V_l, l = 0 \ldots L\). We think of mesh-sizes \(h_l = 2^{-l}\), leading to finite element space dimensions \(N^{dl}\), where \(d\) is the spatial dimension.

We have prolongation matrices

\[ P_l \in {\mathbb R}^{N_l \times N_{l-1}}. \]

If \(v_{l-1}\) is a finite element function in \(V_{l-1}\) represented by the coefficient vector \(\underline v_{l-1}\). Then \(\underline v_l = P_l \underline v_{l-1}\) is the coefficient vector of the same function represented by the basis functions of \(V_l\).

If \(A_l\) and \(A_{l-1}\) are discretization matrices by a Galerkin method, then there holds

\[ A_{l-1} = P_l^T A_l P_l \]

Let \(D_l = \operatorname{diag} A_l\) be the Jacobi preconditioner (or some similar, cheap and local preconditioner).

2-level method

A 2-level preconditioner involving levels \(l-1\) and level \(l\) is

\[ C_{2L}^{-1} = D_l^{-1} + P_l A_{l-1}^{-1} P_l^T \]

By the ASM - Lemma we have the tools to analyze the this preconditioner has optimal condition number. However, a direct inversion of the matrix \(A_{l-1}\) is up to constant factor as expensive as the inversion of \(A_l\). The multilevel method is to replace the inverses recursively by a multilevel preconditioner:

\[\begin{eqnarray*} C_{ML,0}^{-1} & = & A_0^{-1} \\ C_{ML,l}^{-1} & = & D_l^{-1} + P_l C_{ML,l-1}^{-1} P_l^T \qquad \text{for} \; l = 1, \ldots L \end{eqnarray*}\]

from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

fes = H1(mesh,order=1, dirichlet=".*", autoupdate=True)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)

class MLPreconditioner(BaseMatrix):
    def __init__ (self, fes, level, mat, coarsepre):
        super().__init__()
        self.fes = fes
        self.level = level
        self.mat = mat
        self.coarsepre = coarsepre
        if level > 0:
            self.localpre = mat.CreateSmoother(fes.FreeDofs())
        else:
            self.localpre = mat.Inverse(fes.FreeDofs())
        
    def Mult (self, x, y):
        if self.level == 0:
            y.data = self.localpre * x
            return
        hx = x.CreateVector(copy=True)
        self.fes.Prolongation().Restrict(self.level, hx)
        cdofs = self.fes.Prolongation().LevelDofs(self.level-1)
        y[cdofs] = self.coarsepre * hx[cdofs] 
        self.fes.Prolongation().Prolongate(self.level, y)
        y += self.localpre * x

    def Shape (self):
        return self.localpre.shape
    def CreateVector (self, col):
        return self.localpre.CreateVector(col)

With operator notation we can define the multilevel preconditioner also as follows:

def MLPreconditioner2(fes, level, mat, coarsepre):
    prol = fes.Prolongation().Operator(level)
    localpre = mat.CreateSmoother(fes.FreeDofs())
    return localpre + prol @ coarsepre @ prol.T

a.Assemble()
pre = a.mat.Inverse(fes.FreeDofs())

for l in range(8):
    mesh.Refine()
    print ("ndof = ", fes.ndof)
    a.Assemble()
    pre = MLPreconditioner(fes,l+1, a.mat, pre)
    
    lam = EigenValues_Preconditioner(a.mat, pre)
    print ("lammin, lammax=", lam[0], lam[-1], \
            "kappa=", lam[-1]/lam[0])

ndof =  61
lammin, lammax= 0.5225660928343492 1.9634170326035565 kappa= 3.7572606786524703
ndof =  217
lammin, lammax= 0.4164879047972152 3.2550708572760323 kappa= 7.815523139527669
ndof =  817
lammin, lammax= 0.3798934475818364 4.308450122441681 kappa= 11.341206724850283
ndof =  3169
lammin, lammax= 0.3560043386452684 5.15017289227362 kappa= 14.466601479835841
ndof =  12481
lammin, lammax= 0.31993781940106436 5.825268149129326 kappa= 18.20750094513505
ndof =  49537
lammin, lammax= 0.31667701963456707 6.37087318741975 kappa= 20.11788918176472

ndof =  197377

lammin, lammax= 0.31820303480889367 6.81610745286266 kappa= 21.42062364979101

ndof =  787969

lammin, lammax= 0.3219928450515168 7.183205809357794 kappa= 22.308588279992765

58.2. Multigrid Preconditioning

The multigrid preconditioner combines the local preconditioner and the coarse grid step sequentially. The multigrid preconditioning action is define as follows:

Multigrid \(C_{MG,l}^{-1} : d \mapsto w\)

• if l = 0, set \(w = A_0^{-1} d\) and return

• w = 0

• presmoothing, \(m_l\) steps:
  $\(\hat w = w + D_{pre}^{-1} (d - A w)\)$

• coasre grid correction: $\( \hat w = w + P_l C_{MG,l-1}^{-1} P_l^{T} (d - A w) \)$

• postsmoothing, \(m_l\) steps:
  $\(\hat w = w + D_{post}^{-1} (d - A w)\)$

If the preconditioners \(D_{pre}\) and \(D_{post}\) from the pre-smoothing and post-smoothing iterations are transposed to each other, the overall preconditioner is symmetric. If the pre- and post-smoothing iterations are non-expansive, the overall preconditioner is positive definite. This can be obtained by applying forward Gauss-Seidel for pre-smoothing, and backward Gauss-Seidel for post-smoothing.

from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

fes = H1(mesh,order=1, dirichlet=".*")
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)

class MGPreconditioner(BaseMatrix):
    def __init__ (self, fes, level, mat, coarsepre):
        super().__init__()
        self.fes = fes
        self.level = level
        self.mat = mat
        self.coarsepre = coarsepre
        if level > 0:
            self.localpre = mat.CreateSmoother(fes.FreeDofs())
        else:
            self.localpre = mat.Inverse(fes.FreeDofs())
        
    def Mult (self, d, w):
        # print ("level = ", self.level)
        if self.level == 0:
            w.data = self.localpre * d
            return
        
        prol = self.fes.Prolongation().Operator(self.level)

        w[:] = 0
        self.localpre.Smooth(w,d)
        res  = d - self.mat * w
        w += prol @ self.coarsepre @ prol.T * res
        self.localpre.SmoothBack(w,d)


    def Shape (self):
        return self.localpre.shape
    def CreateVector (self, col):
        return self.localpre.CreateVector(col)

a.Assemble()
pre = MGPreconditioner(fes, 0, a.mat, None)

for l in range(5):
    mesh.Refine()
    fes.Update()
    print ("ndof = ", fes.ndof)
    a.Assemble()
    pre = MGPreconditioner(fes,l+1, a.mat, pre)
    
    lam = EigenValues_Preconditioner(a.mat, pre)
    print ("lammin, lammax=", lam[0], lam[-1], \
            "kappa=", lam[-1]/lam[0])

ndof =  61
lammin, lammax= 0.7651544666507635 0.9987999466057931 kappa= 1.305357271163486
ndof =  217
lammin, lammax= 0.5892049430232533 0.9967208414315474 kappa= 1.691636930805944
ndof =  817
lammin, lammax= 0.5497371221971643 0.99706677703541 kappa= 1.8137155683617998
ndof =  3169
lammin, lammax= 0.4814120552195523 0.996951398840777 kappa= 2.070889974672754
ndof =  12481
lammin, lammax= 0.45725397386741373 0.9961888045847007 kappa= 2.1786334543121053

f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)
from ngsolve.krylovspace import CGSolver
inv = CGSolver(mat=a.mat, pre=pre, printrates=True)
gfu.vec.data = inv * f.vec

CG iteration 1, residual = 0.17521646837299754     
CG iteration 2, residual = 0.009411871240583954     
CG iteration 3, residual = 0.0010818082293434137     
CG iteration 4, residual = 0.00017741243435071406     
CG iteration 5, residual = 3.1462473928099643e-05     
CG iteration 6, residual = 5.201427642308935e-06     
CG iteration 7, residual = 7.973354406342009e-07     
CG iteration 8, residual = 1.6349639615836007e-07     
CG iteration 9, residual = 3.6495382905213284e-08     
CG iteration 10, residual = 6.608662509689949e-09     
CG iteration 11, residual = 1.257759608145012e-09     
CG iteration 12, residual = 2.4689978076758385e-10     
CG iteration 13, residual = 4.949419108376888e-11     
CG iteration 14, residual = 9.056463652131037e-12     
CG iteration 15, residual = 1.7024528929216965e-12     
CG iteration 16, residual = 3.143292625080651e-13     
CG iteration 17, residual = 5.494419632353245e-14     

from ngsolve.webgui import Draw
if fes.ndof < 20000:
    Draw(gfu, order=1)

58.3. Projection matrices from the finest level

It is often not feasible to assemble matrices on the coarse level, for example when solving non-linear problems. Then it is useful to calculate coarse grid matrices from the matrix on the finest level using the Galerkin property

\[ A_{l-1} = P_{l}^T A_l P_l \]

[ Requires NGSolve updates from 14.04.2021, still tricky: freedofs on coarser levels]

from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

fes = H1(mesh,order=1, dirichlet=".*", autoupdate=True)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)

for l in range(5):
     mesh.Refine()
    
a.Assemble();

class ProjectedMG(BaseMatrix):
    def __init__ (self, fes, mat, level):
        super(ProjectedMG, self).__init__()
        self.fes = fes
        self.level = level
        self.mat = mat
        if level > 0:
            self.prol = fes.Prolongation().CreateMatrix(level)
            self.rest = self.prol.CreateTranspose()
            coarsemat = self.rest @ mat @ self.prol # multiply matrices
            self.localpre = mat.CreateSmoother(fes.FreeDofs())
                
            self.coarsepre = ProjectedMG(fes, coarsemat, level-1)
        else:
            self.localpre = mat.Inverse(fes.FreeDofs())
        
    def Mult (self, d, w):
        if self.level == 0:
            w.data = self.localpre * d
            return
        w[:] = 0
        self.localpre.Smooth(w,d)
        res = d - self.mat * w
        w += self.prol @ self.coarsepre @ self.rest * res
        self.localpre.SmoothBack(w,d)

        
    def Shape (self):
        return self.localpre.shape
    def CreateVector (self, col):
        return self.localpre.CreateVector(col)

# print (fes.mesh.levels)
pre = ProjectedMG(fes, a.mat, fes.mesh.levels-1)

print ("ndof: ", fes.ndof)
lam = EigenValues_Preconditioner(a.mat, pre)
print ("lammin, lammax=", lam[0], lam[-1], \
            "kappa=", lam[-1]/lam[0])

ndof:  12481
lammin, lammax= 0.4613575315481071 0.9962490853348935 kappa= 2.1593861966270116

f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)
from ngsolve.krylovspace import CGSolver
inv = CGSolver(mat=a.mat, pre=pre, printrates=True)
gfu.vec.data = inv * f.vec

from ngsolve.webgui import Draw
if fes.ndof < 20000:
    Draw(gfu, order=1)

CG iteration 1, residual = 0.17521646837299754     
CG iteration 2, residual = 0.009411871240584117     
CG iteration 3, residual = 0.0010818082293433905     
CG iteration 4, residual = 0.0001774124343507163     
CG iteration 5, residual = 3.1462473928099454e-05     
CG iteration 6, residual = 5.201427642308891e-06     
CG iteration 7, residual = 7.973354406341718e-07     
CG iteration 8, residual = 1.634963961583516e-07     
CG iteration 9, residual = 3.6495382905213066e-08     
CG iteration 10, residual = 6.608662509690066e-09     
CG iteration 11, residual = 1.2577596081449808e-09     
CG iteration 12, residual = 2.468997807675683e-10     
CG iteration 13, residual = 4.949419108376651e-11     
CG iteration 14, residual = 9.05646365213114e-12     
CG iteration 15, residual = 1.7024528929218712e-12     
CG iteration 16, residual = 3.143292625081416e-13     
CG iteration 17, residual = 5.494419632357074e-14